Sensitivity Testing
By: Claire Gonzales
2024-08-10
Read and wrangle in data
output_main <- read.csv(here("data","bivalve_model","LookUpTable_201416.csv")) %>%
dplyr::select("site", "long", "lat", "yield")
output_s1 <- read.csv(here("data","bivalve_model","LookUpTable_s1.csv"))
output_s2 <- read.csv(here("data","bivalve_model","LookUpTable_s2.csv"))
output_s3 <- read.csv(here("data","bivalve_model","LookUpTable_s3.csv"))s1 <- output_s1 %>%
mutate(s1_yield = yield) %>%
dplyr::select("site", "long", "lat", "s1_yield")
s2 <- output_s2 %>%
mutate(s2_yield = yield) %>%
dplyr::select("site", "long", "lat", "s2_yield")
s3 <- output_s3 %>%
mutate(s3_yield = yield) %>%
dplyr::select("site", "long", "lat", "s3_yield")Merging the 3 frames. Format will be wide so then we will make it longer.
wide1 <- merge(output_main, s1, by = c("site", "long", "lat"), all = TRUE)
wide2 <- merge(s2, s3, by = c("site", "long", "lat"), all = TRUE)
wide <- merge(wide1, wide2, by = c("site", "long", "lat"), all = TRUE)alldata <- wide %>%
pivot_longer(yield:s3_yield,
names_to = "test",
values_to = "yield")Visualizations
##boxplot and violin plot
ggplot(data = alldata, aes(x = test, y = yield, fill = test)) +
geom_boxplot() +
#geom_jitter(color="black", size=0.4, alpha=0.9) +
theme(
legend.position="none",
plot.title = element_text(size=11)
) +
theme_minimal()
ggplot(data = alldata, aes(x = test, y = yield, fill = test)) +
geom_violin() +
#geom_jitter(color="black", size=0.4, alpha=0.9) +
theme(
legend.position="none",
plot.title = element_text(size=11)
) +
theme_minimal()
Percent change
wide_percents <- wide %>%
mutate(s1_percchange = (((s1_yield-yield)/yield)*100),
s2_percchange = (((s2_yield-yield)/yield)*100),
s3_percchange = (((s3_yield-yield)/yield)*100),
)
summary_percents <- wide_percents %>%
#group_by(site) %>%
summarize(mean_s1percchhange = mean(s1_percchange),
mean_s2percchhange = mean(s2_percchange),
mean_s3percchhange = mean(s3_percchange)
)ANOVA
summary(alldata)## site long lat test
## Min. : 6 Min. :-124.6 Min. :32.35 Length:18396
## 1st Qu.:1240 1st Qu.:-122.8 1st Qu.:33.26 Class :character
## Median :2415 Median :-119.7 Median :33.93 Mode :character
## Mean :2476 Mean :-120.5 Mean :35.33
## 3rd Qu.:3757 3rd Qu.:-118.7 3rd Qu.:37.67
## Max. :5000 Max. :-117.1 Max. :42.01
## yield
## Min. : 871687
## 1st Qu.: 2836668
## Median : 3634054
## Mean : 4433778
## 3rd Qu.: 5955865
## Max. :11570374aov <- aov(yield ~ test,
data = alldata)
summary(aov)## Df Sum Sq Mean Sq F value Pr(>F)
## test 3 2.520e+16 8.400e+15 2831 <2e-16 ***
## Residuals 18392 5.458e+16 2.968e+12
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1hist(aov$residuals)
ggplot(data = alldata, aes(sample = yield)) +
geom_qq() +
facet_wrap( ~ test) +
theme_minimal()
Visualization data transformations
This section deals with the production data and transformed production data from each of the 4 sub-models. We want to compare the transformed data to the theoretical transformation (line resulting from the equation we used to transform the data - how well do the data points sit on this line?).
I also want to use this opportunity to look at the production data (pre-transformation) in some non-spatial ways.
I will read in the rasters and then extract the value in each cell.
Read in data using terra
library(terra)
library(sf)Wave energy
wave <- terra::rast(here("data", "tiff", "wave_linfunc_v1.tif"))
wave_raw <- terra::rast(here("data", "tiff", "wave_notransform.tif"))
terra::plot(wave, main = "Wave Energy (transformed linearly)")
terra::plot(wave_raw, main = "Wave Energy (no transformation)")
View head and tail of data frame
wave_df <- terra::as.data.frame(wave, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)
head(wave_df)## x y Band_1
## 4884 -124.5397 41.97784 0.8077396
## 4886 -124.4588 41.97784 0.7711606
## 4887 -124.4184 41.97784 0.7387354
## 4888 -124.3780 41.97784 0.6629785
## 4889 -124.3376 41.97784 0.6135752
## 4890 -124.2971 41.97784 0.5418847tail(wave_df)## x y Band_1
## 51040 -119.0416 32.43699 0.6873534
## 51041 -119.0012 32.43699 0.6822757
## 51235 -119.0416 32.39656 0.6925157
## 51236 -119.0012 32.39656 0.6886186
## 51237 -118.9607 32.39656 0.6806911
## 52815 -118.2330 32.07314 0.5573189ggplot(data = wave_df, aes(x = Band_1)) +
geom_histogram() +
theme_minimal()
###################################################
wave_raw_df <- terra::as.data.frame(wave_raw, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)
head(wave_raw_df)## x y wave_notransform
## 4884 -124.5397 41.96992 38814.79
## 4886 -124.4588 41.96992 36715.67
## 4887 -124.4184 41.96992 34854.92
## 4888 -124.3780 41.96992 30507.54
## 4889 -124.3376 41.96992 27672.48
## 4890 -124.2971 41.96992 23558.46tail(wave_raw_df)## x y wave_notransform
## 51040 -119.0416 32.43665 31906.31
## 51041 -119.0012 32.43665 31614.92
## 51235 -119.0416 32.39625 32202.55
## 51236 -119.0012 32.39625 31978.92
## 51237 -118.9607 32.39625 31523.99
## 52815 -118.2330 32.07309 24444.17ggplot(data = wave_raw_df, aes(x = wave_notransform)) +
geom_histogram() +
theme_minimal()
# Creating a scatterplot of data with the transformation line on it
graph_wa <- ggplot(data = wave_df, aes(x = , y =)) +
geom_point() +
theme_minimal()Merging the pre-transformed and transformed data frames
####### YIKES THIS DOESNT WORK ################
# wave_all <- merge(wave_raw_df, wave_df, by = c("x", "y"), all = FALSE)Wind energy
wind <- terra::rast(here("data", "tiff", "wind_linfunc_v_ProjectRas.tif"))
wind_raw <- terra::rast(here("data", "tiff", "wind_notransform.tif"))
terra::plot(wind, main = "Wind Energy (transformed linearly)")
terra::plot(wind_raw, main = "Wind Energy (no transformation)")
View head and tail of data frame
wind_df <- terra::as.data.frame(wind, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)
head(wind_df)## x y Band_1
## 1083981 -124.3807 42.06418 0.7324007
## 1083984 -124.3708 42.06418 0.7286220
## 1083985 -124.3675 42.06418 0.7284269
## 1083988 -124.3577 42.06418 0.7237991
## 1083989 -124.3544 42.06418 0.7185299
## 1083991 -124.3478 42.06418 0.7152024tail(wind_df)## x y Band_1
## 10049327 -119.0085 32.16150 0.4130057
## 10052294 -119.0381 32.15821 0.4195321
## 10052295 -119.0348 32.15821 0.4183595
## 10052296 -119.0315 32.15821 0.4169410
## 10052297 -119.0282 32.15821 0.4160608
## 10052298 -119.0250 32.15821 0.4160352ggplot(data = wind_df, aes(x = Band_1)) +
geom_histogram() +
theme_minimal()
###################################################
wind_raw_df <- terra::as.data.frame(wind_raw, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)
head(wind_raw_df)## x y Band_1
## 990023 -124.3441 42.06536 23472.07
## 990026 -124.3342 42.06536 23263.02
## 990027 -124.3309 42.06536 23065.66
## 990029 -124.3243 42.06536 23032.67
## 990030 -124.3210 42.06536 22990.36
## 990031 -124.3177 42.06536 22983.74tail(wind_raw_df)## x y Band_1
## 9913185 -119.0245 32.16039 15022.49
## 9913186 -119.0212 32.16039 15014.09
## 9913187 -119.0179 32.16039 15005.87
## 9913188 -119.0146 32.16039 14989.09
## 9916144 -119.0343 32.15710 15081.72
## 9916145 -119.0310 32.15710 15048.69ggplot(data = wind_raw_df, aes(x = Band_1)) +
geom_histogram() +
theme_minimal()
Merging
# wind_all <- merge(wind_raw_df, wind_df, by = c("x", "y"), all = FALSE)Aquaculture
aqua <- terra::rast(here("data", "tiff", "aqua_sfunc_v2.tif")) #???????
aqua_raw <- terra::rast(here("data", "tiff", "aqua_notransform.tif"))
terra::plot(aqua, main = "Aquaculture Production (S-transformation)")
terra::plot(aqua_raw, main = "Aquaculture Production (no transformation)")
View head and tail of data frame
aqua_df <- terra::as.data.frame(aqua, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)
head(aqua_df)## x y Band_1
## 8693 -124.4787 42.00743 0.5604585
## 8694 -124.4487 42.00743 0.6017970
## 8695 -124.4187 42.00743 0.6330283
## 8696 -124.3888 42.00743 0.6407387
## 8697 -124.3588 42.00743 0.6946170
## 8698 -124.3288 42.00743 0.7214648tail(aqua_df)## x y Band_1
## 93034 -119.0533 32.41554 4.251122e-04
## 93035 -119.0233 32.41554 8.228081e-05
## 93036 -118.9933 32.41554 9.777889e-06
## 93297 -119.0533 32.38557 1.366302e-04
## 93299 -118.9933 32.38557 2.556847e-06
## 93300 -118.9633 32.38557 1.362262e-06ggplot(data = aqua_df, aes(x = Band_1)) +
geom_histogram() +
theme_minimal()
###################################################
aqua_raw_df <- terra::as.data.frame(aqua_raw, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)
head(aqua_raw_df)## x y aqua_notransform
## 8693 -124.4787 41.99691 7521588
## 8694 -124.4487 41.99691 7716453
## 8695 -124.4187 41.99691 7870482
## 8696 -124.3888 41.99691 7909508
## 8697 -124.3588 41.99691 8194795
## 8698 -124.3288 41.99691 8346378tail(aqua_raw_df)## x y aqua_notransform
## 93034 -119.0533 32.41518 3064938
## 93035 -119.0233 32.41518 2994500
## 93036 -118.9933 32.41518 2958244
## 93297 -119.0533 32.38523 3010470
## 93299 -118.9933 32.38523 2948923
## 93300 -118.9633 32.38523 2946289ggplot(data = aqua_raw_df, aes(x = aqua_notransform)) +
geom_histogram() +
theme_minimal()
Merging
aqua_all <- merge(aqua_raw_df, aqua_df, by = c("x", "y"), all = FALSE)Fisheries
SOMETHNG IS STILL WRONG HERE. I think maybe there might be an issue with the pre-transformation raster still. Transformed data has 1,640,319 observations and pre-transformation only has 15,541, which is weird.
fishy <- terra::rast(here("data", "tiff", "catchdata_outp_ProjectRas.tif"))
fishy_raw <- terra::rast(here("data", "tiff", "catchdata_notransform.tif"))
terra::plot(fishy, main = "Fisheries Production (Z-transformation)")
terra::plot(fishy_raw, main = "Fisheries Production (no transformation)")
View head and tail of data frame
fishy_df <- terra::as.data.frame(fishy, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)
head(fishy_df)## x y catchdata_outp_ProjectRas
## 940787 -125.5002 41.99897 0.6294906
## 940788 -125.4969 41.99897 0.6294906
## 940789 -125.4936 41.99897 0.6294906
## 940790 -125.4903 41.99897 0.6294906
## 940791 -125.4869 41.99897 0.6294906
## 940792 -125.4836 41.99897 0.6294906tail(fishy_df)## x y catchdata_outp_ProjectRas
## 9273852 -117.1849 32.33468 0.8797805
## 9273853 -117.1816 32.33468 0.8797805
## 9273854 -117.1782 32.33468 0.8797805
## 9273855 -117.1749 32.33468 0.8797805
## 9273856 -117.1716 32.33468 0.8797805
## 9273857 -117.1683 32.33468 0.8797805ggplot(data = fishy_df, aes(x = catchdata_outp_ProjectRas)) +
geom_histogram() +
theme_minimal()
###################################################
fishy_raw_df <- terra::as.data.frame(fishy_raw, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)
head(fishy_raw_df)## x y catchdata_notransform
## 7 -125.479 41.99396 25653.79
## 8 -125.445 41.99396 25653.79
## 9 -125.411 41.99396 25653.79
## 10 -125.377 41.99396 25653.79
## 11 -125.343 41.99396 25653.79
## 12 -125.309 41.99396 25653.79tail(fishy_raw_df)## x y catchdata_notransform
## 72054 -118.849 32.33796 211535.130
## 72099 -117.319 32.33796 2010.012
## 72100 -117.285 32.33796 2010.012
## 72101 -117.251 32.33796 2010.012
## 72102 -117.217 32.33796 2010.012
## 72103 -117.183 32.33796 2010.012ggplot(data = fishy_raw_df, aes(x = catchdata_notransform)) +
geom_histogram() +
theme_minimal()
# xlim(c(-1, 500000))fishy_all <- merge(fishy_raw_df, fishy_df, by = c("x", "y"), all = FALSE)Co-location Suitability Model
colo <- terra::rast(here("data", "tiff", "CellSta_colocation_off.tif"))
terra::plot(colo, main = "Co-location Suitability Score")
colo_df <- terra::as.data.frame(colo, row.names = NULL, optional = NULL, xy = TRUE, cells=FALSE, time=FALSE, na.rm=TRUE, wide=TRUE)
head(colo_df)## x y Band_1
## 490 -125.4974 41.98627 0.6294906
## 491 -125.4570 41.98627 0.7640004
## 492 -125.4165 41.98627 0.7640711
## 493 -125.3761 41.98627 0.7652843
## 494 -125.3357 41.98627 0.7650780
## 495 -125.2952 41.98627 0.7662386tail(colo_df)## x y Band_1
## 56103 -119.0694 32.36456 0.3808336
## 56104 -119.0290 32.36456 0.3760225
## 56105 -118.9886 32.36456 0.3719020
## 56106 -118.9482 32.36456 0.3670245
## 56107 -118.9077 32.36456 0.3618057
## 56108 -118.8673 32.36456 0.3588376ggplot(data = colo_df, aes(x = Band_1)) +
geom_histogram() +
theme_minimal()
Scatter plot showing the co-location suitability score related to latitude
# library(hrbrthemes) #having issues downloading some of these dependencies
ggplot(data = colo_df, aes(x = y, y = Band_1)) +
geom_point(size = 1) +
geom_smooth(method="lm") +
theme_minimal()
ggplot(data = colo_df, aes(x = x, y = Band_1)) +
geom_point(size = 1) +
geom_smooth(method="lm") +
theme_minimal()
Linear models to correspond each graph above
long_lm <- lm(Band_1 ~ y, data = colo_df)
lat_lm <- lm(Band_1 ~ x, data = colo_df)
# library(broom)
# broom::tidy(long_lm)
# broom::tidy(lat_lm)
# glance(long_lm)
# glance(lat_lm)
summary(long_lm)##
## Call:
## lm(formula = Band_1 ~ y, data = colo_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.50800 -0.10964 0.02131 0.12494 0.49704
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.4783852 0.0215951 -22.15 <2e-16 ***
## y 0.0300670 0.0005907 50.90 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1767 on 10924 degrees of freedom
## Multiple R-squared: 0.1917, Adjusted R-squared: 0.1916
## F-statistic: 2591 on 1 and 10924 DF, p-value: < 2.2e-16summary(lat_lm)##
## Call:
## lm(formula = Band_1 ~ x, data = colo_df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.51035 -0.10076 0.00483 0.10449 0.51011
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -5.9849297 0.0798722 -74.93 <2e-16 ***
## x -0.0540747 0.0006541 -82.68 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1541 on 10924 degrees of freedom
## Multiple R-squared: 0.3849, Adjusted R-squared: 0.3848
## F-statistic: 6835 on 1 and 10924 DF, p-value: < 2.2e-16Latitude seems to have the better correlation and higher R value
plot(lat_lm)
Visualizations
# wave_graph <- ggplot(data = wave_all, aes(x = Band_1, y = wave_notransform)) +
# geom_point() +
# theme_minimal()